MULTINEST is a Bayesian posterior sampler that has two distinct advantages over traditional MCMC:
To run the MULTINEST sampler in 3ML, one must have the foloowing software installed:
MULTINEST can be run in a single instance, but it can be incredibly slow. Luckily, it can be built with MPI support enabling it to be run on a multicore workstation or cluster very effeciently.
There are multiple ways to invoke the parallel run of MULTINEST in 3ML: e.g., one can write a python script with all operations and invoke:
$> mpiexec -n <N> python my3MLscript.py
However, it is nice to be able to stay in the Jupyter environment with ipyparallel which allow the user to easily switch bewteen single instance, desktop cores, and cluster environment all with the same code.
The user is expected to have and MPI distribution installed (open-mpi, mpich) and have compiled MULTINEST against the MPI library. Additionally, the user should setup and ipyparallel profile. Instructions can be found here: http://ipython.readthedocs.io/en/2.x/parallel/parallel_mpi.html
Details for luanching ipcluster on a distributed cluster are not covered here, but everything is the same otherwise.
In the directory that you want to run 3ML in the Jupyter notebook launch and ipcontroller:
$> ipcontroller start --profile=mpi --ip='*'
Next, launch MPI with the desired number of engines:
$> mpiexec -n <N> ipengine --mpi=mpi4py --profile=mpi
Now, the user can head to the notebook and begin!
In [1]:
from ipyparallel import Client
rc = Client(profile='mpi')
# Grab a view
view = rc[:]
# Activate parallel cell magics
view.activate()
Import 3ML and astromodels to the workers
In [2]:
with view.sync_imports():
import threeML
import astromodels
Now we set up the analysis in the normal way except the following two caveats:
In [3]:
%%px
# Make GBM detector objects
src_selection = "0.-10."
nai0 = threeML.FermiGBM_TTE_Like('NAI0',
"glg_tte_n0_bn080916009_v01.fit",
"-10-0,100-200", # background selection
src_selection, # source interval
rspfile="glg_cspec_n0_bn080916009_v07.rsp")
nai3 = threeML.FermiGBM_TTE_Like('NAI3',"glg_tte_n3_bn080916009_v01.fit",
"-10-0,100-200",
src_selection,
rspfile="glg_cspec_n3_bn080916009_v07.rsp")
nai4 = threeML.FermiGBM_TTE_Like('NAI4',"glg_tte_n4_bn080916009_v01.fit",
"-10-0,100-200",
src_selection,
rspfile="glg_cspec_n4_bn080916009_v07.rsp")
bgo0 = threeML.FermiGBM_TTE_Like('BGO0',"glg_tte_b0_bn080916009_v01.fit",
"-10-0,100-200",
src_selection,
rspfile="glg_cspec_b0_bn080916009_v07.rsp")
# Select measurements
nai0.set_active_measurements("10.0-30.0", "40.0-950.0")
nai3.set_active_measurements("10.0-30.0", "40.0-950.0")
nai4.set_active_measurements("10.0-30.0", "40.0-950.0")
bgo0.set_active_measurements("250-43000")
# Set up 3ML likelihood object
triggerName = 'bn080916009'
ra = 121.8
dec = -61.3
data_list = threeML.DataList( nai0,nai3,nai4,bgo0 )
band = astromodels.Band()
GRB = threeML.PointSource( triggerName, ra, dec, spectral_shape=band )
model = threeML.Model( GRB )
# Set up Bayesian details
bayes = threeML.BayesianAnalysis(model, data_list)
band.K.prior = astromodels.Log_uniform_prior(lower_bound=1E-2, upper_bound=5)
band.xp.prior = astromodels.Log_uniform_prior(lower_bound=1E2, upper_bound=2E3)
band.alpha.prior = astromodels.Uniform_prior(lower_bound=-1.5,upper_bound=0.)
band.beta.prior = astromodels.Uniform_prior(lower_bound=-3.,upper_bound=-1.5)
Out[3]:
Finally we call MULTINEST. If all is set up properly, MULTINEST will gather the distributed objects and quickly sample the posterior
In [4]:
%px samples = bayes.sample_multinest(n_live_points=400,resume=False)
Out[4]:
Now we need to bring the BayesianAnalysis object back home. Unfortunately, not all objects can be brought back. So you must save figures to the workers. Future implementations of 3ML will allow for saving of the results to a dedicated format which can then be viewed on the local machine. More soon!
In [17]:
# Execute commands that allow for saving figures
# grabbing samples, etc
%%px --targets ::1
samples = bayes.raw_samples()
f=bayes.get_credible_intervals()
bayes.corner_plot(plot_contours=True, plot_density=False)
Out[17]:
In [22]:
# Bring the raw samples local
raw_samples=view['samples'][0]
In [23]:
raw_samples['bn080916009.spectrum.main.Band.K']
Out[23]:
In [ ]: